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ABSTRACT 

This  paper  considers  the  deterministic  and  stochastic  modeling  of  2-D 
systems  described  by  their  input/output  data.  In  the  deterministic  case, 
the  modeling  problem  is  formulated  as  a  2-D  Pade  approximation  problem.  By 
studying  several  possible  geometries  of  approximation,  we  obtain  several 
sets  of  recursions  of  the  2-D  rational  approximants .  These  results  exploit 
the  properties  of  2-D  Hankel  matrices,  and  they  are  used  here  to  characterize 
2-D  rational  transfer  functions.  In  the  stochastic  case,  the  realization 
problem  is  viewed  as  a  2-D  prediction  problem.  This  problem  is  solved  re¬ 
cursively  by  generalizing  to  the  2-D  case  an  algorithm  due  to  Levinson  in 
the  1-D  case.  The  predictors  obtained  by  this  algorithm  are  then  showed  to 
converge  to  the.  2-D  spectral  factors  of  the  output  spectrum. 
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I.  Introduction 


In  this  paper,  the  modeling  of  2-D  linear  and  shift-invariant  systems  given 
by  their  input/output  data  will  be  considered.  This  problem  will  be  studied  both 
in  the  deterministic  case,  i.e.  when  the  point  spread  function  of  the  system  is 
known,  and  in  the  stochastic  case  when  only  input  and  output  covariances  are  given. 
For  convenience,  we  shall  restrict  our  attention  to  the  case  of  single  input- 
single  output  systems. 

In  the  deterministic  case,  the  realization  problem  will  be  formulated  as  a 
2-D  Pade  approximation  problem.  Thus,  given  a  point  spread  function 

g(z,ui)  =  g  z"V3  (1) 

i  ,j 


we  shall  approximate  g  by  rational  functions 


«*•“>  *  SfeS- 


(2) 


where  a  and  b  are  2-D  polynomials.  Such  a  problem  has  been  considered  by  Chisholm  [1 
Hughes  Jones  and  Makinson  [2],  Graves -Morris,  Hughes  Jones  and  Makinson  [3]  among  oth 

but  the  results  discussed  here  will  be  somewhat  different.  The  main  difference  is 
that  the  realization  algorithms  presented  in  Sections  II  and  III  are  recursive.  . 

These  algorithms  extend  to  the  2-D  case  a  realization  procedure  originally  intro¬ 
duced  by  Lanczos  [4],  Berlekamp  [5],  Massey  [6]  and  Rissanen  [7]  in  the  1-D  case. 

The  recursions  that  we  obtain  exploit  the  properties  of  Hankel  block  Hankel  matrices; 
These  properties  can  be  used  to  relate  the  recursions  obtained  for  g  to  those  de¬ 
rived  by  Jackson  [8]  for  2-D  orthogonal  polynomials  on  the  hyper  real  line. 

In  Section  IV  we  study  the  relation  existing  between  the  2-D  partial  realiza¬ 
tion  problem  considered  here  and  the  complete  realization  problem  studied  by  Ho 
and  Kalman  [9],  Silverman  [10]  and  Youla  and  Tissi  [11]  in  the  1-D  case,  and  by  1 
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Fliess  [12],  Fornasini  [13],  Clerget  [14]  and  Kao  and  Chen  [15]  in  the  2-D  case. 

By  doing  so,  we  obtain  a  set  of  simple  conditions  on  the  Markov  parameters  of  g 
that  can  be  used  to  verify  whether  the  partial  realization  §  =  b/a  is  also  a  com¬ 
plete  realization,  i.e.  whether  g  =  g. 

The  stochastic  realization  problem  is  considered  in  Section  V  and  is  formu¬ 
lated  as  an  autoregressive  modeling  problem.  Thus,  if  the  system  considered  is 
linear  and  shift  invariant,  and  if  it  is  driven  by  a  2-D  white  noise  process  u(i,j), 
the  output  process  will  be  modeled  as 

y(i,j)  +  y  a  (k,£)  y(i-k,j-£)  =  u(i,j)  (3) 

I-(0,0) 

where  I  is  a  causal  asymetric  half  plane  set  of  the  type  considered  in  [ 16] - [ 18] . 

The  coefficients  aj(k,£)  of  the  filter  (3)  will  be  selected  so  that 

y(ij)  3  ~y  a  (k,£)  y(i-k, j-£)  (4) 

I-CO.O) 

is  the  linear  least-squares  estimate  of  y(i,j)  given  observations  over  the  set 
=  { (i-k, j-£)  :  (k,£)  el-  (0,0)}. 

To  solve  this  linear  least-squares  estimation  problem,  we  shall  use  a  2-D  version 
of  an  algorithm  originally  introduced  by  Levinson  [19]  in  the  1-D  case.  The  re¬ 
cursions  that  we  obtain  for  aj(.,.)  were  first  described  in  [20],  [21],  and  they 
differ  from  those  presented  by  Justice  [22]  by  the  fact  that  the  orders  n  and  m  of 
a^  in  z  and  w  can  be  increased  separately  (Justice's  recursions  were  requiring  that 
either  n  or  m  be  fixed  a  priori).  These  recursions  present  also  some  similarity 
with  the  recursions  derived  by  Genin  and  Kamp  [23]  for  2-D  orthogonal  polynomials 
on  the  unit  hypercircle  (see  [17] -[18]  and  [24] -[27]  for  more  details). 


In  addition,  we  relate  the  stochastic  modeling  problem  to  the  spectral 
factorization  results  of  Helson  and  Lowdenslager  [28],  Pistor  [29],  Ekstrom  and 
Woods  [16],  Murray  [30]  and  Genin  and  Kamp  [18].  It  is  shown  in  this  context 
that  when  the  domain  I  becomes  infinite,  aj(z,o))  converges  to  the  spectral  factor 
of  the  output  spectrum  r(z,w).  Finally,  the  Section  VI  contains  some  observa¬ 
tions  on  some  open  problems  and  on  possible  extensions  of  the  results  discussed 


II.  The  2-D  Pade  Approximation  Problem 


Throughout  the  following  sections,  it  will  be  assumed  that  the  2-D  transfer 
function  g(z,w)  that  we  want  to  realize  is  a  South-West  (SW)  quarter  plane  causal 
filter,  i.e. 

8(2,oi)  =  g„  z'V^  (5) 

i>0,  j>0 

so  that  g(z,o))  will  be  modeled  by  a  transfer  function  g  =»  b/a  such  that 

n  m 

a(i,U)  •  J  ^  a. .  ,  a00  =  1  (6a) 

i=0  j=0 


is  monic,  and 

b(z,oi) 


n-i  m- j 
z  0)  J 


(6b) 


There  is  no  loss  of  generality  in  making  the  previous  assumption  since  an  arbitrary 
transfer  function  g(z,oi)  can  always  be  decomposed  into  four  parts  which  are  causal 
in  each  of  the  four  quadrants,  i.e. 

g  =  gSW  +  gNW  +  gSE  +  gNE' 

These  four  parts  can  then  be  approximated  separately. 

In  the  following,  we  shall  consider  two  types  of  Pade  approximants  for  g. 

Definition  Let  §(z,o))  =  b(z,u>)/a(z,u))  be  a  2-D  rational  function.  Then  g  is  said 
to  be  a  rational  PadS  approximant  of  g  in  the  domain  D  c  IN  if 


/ 

/ 

- -  . 

•  t 
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g(z»w)  -  §(z,<d)  =  ^  d  z_iof3  (7) 

(i,j)  £  D 

where  D  is  the  complement  of  D.  Similarly,  g  is  said  to  be  a  modified  Pad£  approxi- 
mant  of  g  over  the  domain  A  C  2^  if 

a(z,w)  g(z,o))  -  b(z,o»)  =  ^  d^z^uf3  (8) 

(i.j)  £  A 

where  A  is  the  complement  of  A. 

The  motivation  for  making  the  distinction  between  the  problems  (7)  and  (8)  is 
that  it  is  not  always  possible  to  convert  a  rational  approximation  problem  into 
a  modified  one  and  vice-versa.  An  example  of  domain  of  modified  approximation 
that  cannot  be  converted  into  an  equivalent  domain  of  rational  approximation  is 
described  in  Figure  1.  In  this  case,  the  parameters  of  afz.u)  and  b(z,co)  satisfy 
2(n+l)(m+l)  linear  constraints,  but  if  we  consider  g  =  b/a,  the  Markov  parameters 
g^  of  g(z,w)  are  matched  by  those  of  g  only  over  the  domain 

D(n,m)  ■  { (i, j)  :  (0,0)  <  (i,j)  <  (n,m)} 

which  contains  (n+l)(m+l)  parameters.  The  modified  approximation  problem  associated 
to  Figure  1  will  be  considered  in  Section  III.  However,  we  shall  consider  also 
several  geometries  of  rational  approximation  that  can  be  converted  easily  into  modi¬ 
fied  ones.  The  Figure  2  describes  several  cases  that  will  be  discussed  below. 

To  justify  the  study  of  both  the  problems  (7)  and  (8) ,  it  is  useful  to  note 
that  when  g  =  b/a  is  rational,  the  rational  approximant  gR  and  the  modified  approxi- 
raant  gw  are  such  that 

M 


2r  *  ■  b'a 

provided  that  the  domains  D  and  A  are  chosen  sufficiently  large.  In  the  remainder 
of  this  section  we  shall  discuss  the  rational  approximation  problem  for  the  cases 
A-C  of  Figure  2. 

Case  A:  This  case  was  the  one  considered  by  Chisholm  [1],  Hughes  Jones  and  Makin- 
son  [2]  and  Graves-Morris,  Hughes  Jones  and  Makinson  [3].  In  this  case,  the  Pade 
approximants  have  several  interesting  properties  (such  as  invariance  under  a  large 
class  of  transformations) ,  but  this  geometry  does  not  permit  the  efficient  exploit¬ 
ation  of  the  shift-invariance  properties  of  the  realization  problem.  This  means 
that  the  Pade  approximants  1[n  n  (z  ,u0  where 

n  *  degz  a  *  deg^  a 

cannot  be  computed  recursively. 

Case  B:  If  the  domain  of  rational  approximation  is  given  by 

D(p,q)  =  {(i,j)  :  (0,0)  <  (i,j)  <  (p,q)}  (9) 

and  if  we  assume  that 

(n  =  degz  a,  m  =  deg^  a)  <  (p,q) 

then  the  rational  approximation  problem  (7)  can  be  transformed  into  an  equivalent 
modified  approximation  problem  where  the  domain  A  of  approximation  is  given  by 

A(n,m;  p,q)  =  {(i,j)  :  (-n,-m)  <  (i,j)  <  (p-n,  q-m)}.  (10) 

This  means  that  if  we  denote  by  the  projection  operator  which  selects  the  co- 


k  « 


-v  i  i i«w ii 


efficients  of  a  Laurent  power  series 
h(z,u))  =^  .  z^oT^ 

1iZ 

which  belong  to  A,  i.e. 

nAh(z,uO  *  ^  hf .  z'Vj  , 

A 

then  a(zfoi)  and  b(z,u)  satisfy  the  equation 

nACn,m;  p,q)  *  b(z'“>  =  0  *  CU) 

Consequently,  the  coefficients  of  a(z,to)  and  b(z,o))  obey  a  set  of  linear 
equations  given  by 

bij  *  I  h-Ki-t  V-  e  (12) 

D(n,m) 


2  gi_k, j-f  ak^*  e  D(P’^  *  D(n»m)  •  CW) 

D(n,m) 


By  scanning  the  set  D(p,q)  row  by  row  and  by  denoting 


aj (n;p)  - 


aoj 

b  . 

03 

aij 

b.(n;p)  = 

b.  . 
ij 

anj. 

b 

ng 

°p-n  | 

1 

0 

p-n 

°  1  3  1  m  (14a) 


/ 
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GjCp)  = 


iJ'- 


P] 


g.  .  g  . 
Sl]  6oj 


(14b) 


where  0  is  a  zero  p-n  *  1  vector,  we  can  rewrite  (12)  and  (13)  in  matrix  form, 
p-n 

Indeed,  if  ®  denotes  the  matrix  Kronecker  product,  and  if  we  define 


a(n,ra;  p.q)  = 


aQ(n;p) 


0  so, 

q-m  p+1 


b(n,m.;  p.q)  = 


bQ(n;p)  \ 


b(n;p) 

m 

0  0  0, 
q-m  p+1 


'(P.q) 


G0(P) 


G.  (p) 


Gq(p) 


\ 


\ 


\ 


\ 


\ 


\ 


x 


X 

"  G^p)  Go'  p) 


(15) 


it  is  easy  to  see  that  (12)  and  (13)  are  equivalent  to 


b(n,ra;  p.q)  =  G(p,q)  a(n,m;  p,q)  . 


(16) 


This  equation  has  several  interesting  features.  One  of  them  is  that  G(p,q)  is 

a  block  lower  triangular  Toeplitz  matrix  whose  blocks  are  themselves  lower  trian- 

2 

gular  Toeplitz.  This  property  will  be  denoted  as  G(p,q)  e  LT  ,  and  it  will  be 


shown  below  that  this  structure  has  several  advantages.  Another  interesting  aspect 


-10- 


of  (16)  is  that  the  coefficients  of  b(z,w)  are  entirely  given  by  those  of  a(.,.), 
so  that  the  main  problem  is  to  compute  a(.,.). 

To  solve  (16),  one  needs  to  match  the  (p+1) (q+1)  elements  of  A(n,m;  p,q)  with 
the  2(n+l)(n+l)  coefficients  of  a(.,.)  and  b  (.,.),  so  that  we  have  to  select  p  and 
q  such  that  2(n+l)(m+l)  (p+1)  (q+1).  One  such  choice  is  given  by 

p  =  2n  +  l,q=»m  (17a) 

or  symmetrically  by 

p  =  n  ,  q  =  2m  +  1  (17b) 

This  choice  corresponds  to  the  geometry  C  of  Figure  2.  Other  choices  are  possible, 
but  they  do  not  seem  to  give  rise  to  recursions  in  n  and  m  for  the  polynomials 
a(.,.)  and  b(.,.).  Another  approach  that  will  be  considered  in  Section  III  is 
when 

p  =  2n  +  l  ,  q  =  2m+l  (18) 

and  when  only  the  parameters  of  the  subset  of  A(n,m;  2n+l,  2m+l)  described  in 
Figure  1  are  matched.  As  we  shall  see,  this  case  is  the  most  closely  related  to 
the  1-D  case  (it  depends  on  the  properties  of  2-D  Hankel  matrices) . 


Case  C:  If  p  and  q  are  selected  as  in  (17a) ,  we  can  denote  by 


&  ( z ,  uj)  = 

5n,m  v  ’  J  a 


n,m 


(z,u) 


n,m 


(z,tu) 


the  Pad6  approximant  that  matches  the  Markov  parameters  g^  in  the  domain  D(2n+l,m) 

Since  a  (z.w)  is  chosen  to  be  monic  (i.e.  its  leading  coefficient  a_  _  (0,0)=1), 
n,m  — . — 

the  number  of  free  coefficients  of  §  (z,<*))  is  only  2(n+l)  (m+l)-l,  which  is  one 

n  ,m 


less  than  the  number  of  parameters  in  D(2n+1,  m)  or  in  A(n,m;  2n+l,  m) .  Therefore, 
it  will  be  assumed  in  the  following  that  the  coefficient  of  z"^n+1^a)m  in  A(n,m); 
2n+l,  m)  is  not  matched.  In  this  case,  (11)  becomes 

IT.,  _  ,  ,  a  (z,oj)  g(z,u>)  =  b  (z,oj)  +5  z  ^n+^ojm  (19) 

A(n,m;  2n+l,  m)  n,m  n,m  ’  1  n,m  v  J 

where  <$n  is  the  residual  corresponding  to  the  mis-match  of  z  ^n+^com. 

Since  A(n,m;  2n+l,  m)  =  -D(n,m)U  2(n,m),  where 


-D(n,m)  =  {(i, j)  :  (-i,-j)  e  D(n,m)} 


2(n,m)  =  { (i, j)  :  (i-l,-j)  e  D(n,m)} 


we  can  decompose  (19)  into  two  parts: 


"-OCn.m)  an,m  (z'“>  2 •“>*  bn,m  tz’“> 


expresses  b  (z,ui)  in  function  of  a  (z,go),  and 
n,m  n,m 

n_,  ,  a  (z,u)  g(z,u>)  =  6  z  Cn+1)0)m 

2(n,*m)  n,m  v  J  ’  n,m 


is  the  equation  satisfied  by  a^  m  (z,u)). 
tions,  one  can  denote  by 


aQ (n,m) 


To  obtain  the  matrix  form  of  these  equa- 


bQ(n,m) 


a(n,m)  *  I  a..(n,m) 


a  (n,m) 
m 


an,n/°’^ 

an,n,CI'j) 

an,m(n'j) 


b(n,m)  =  b.  (n,m) 


b  (n,m) 
m 


b  _(o,j) 

n  ,m 

b  (i» j) 

n,m 

b  (n, j) 
n,m 


a^  (n,m) 


bj  (n,m)  * 


the  matrix  form  of  (21)  is  given  by 


x(n,m)  a(n,m)  =  6(n,m)  .  (26) 

To  compute  a(n,m)  and  b(n,m)  recursively,  we  will  now  exploit  the  structure  of 
the  matrices  G(n,m)  and  T(n,m).  This  gives  the  following  recursions. 


Increase  in  m 

The  lower  triangular  structure  of  T(n,m)  can  be  exploited  by  assuming  that 


in  (22)  the  vectors  a^ (n,m)  do  not  depend  on  m,  i.e. 


a^.(n,m)  =  a.  (n)  for  all  j. 


(27) 


This  implies  that  the  vectors  b^(n,m)  do  not  depend  on  m  either.  Thus,  given 


a(n,m)  and  b(n,m),  to  compute  a(n,m+l)  and  b(n,m+l),  we  need  only  to  find  am+^(n) 


and  bm+1(n).  This  can  be  done  by  direct  substitution,  so  that 

m 


Wn)  =  T01  GO  2  Vl-j(n) 

j=i 

ra+l 

bm+lfn)  =  2  Gj(n)  Vl-jCn)- 

j+0 


(28) 


(29) 


■Since  TQ(n)  is  a  Toeplitz  matrix,  its  inverse  can  be  computed  with  0(n  )  opera¬ 
tions  by  using  the  inversion  algorithm  of  Levinson  [19]  and  Trench  [31],  or  with 
2 


0(n  log  n)  operations  if  we  use  a  more  efficient  version  of  this  algorithm  based 
on  doubling  ideas  (cf.  Gustavson  and  Yun  [32],  Morf  [33]  and  Bitmead  and  Anderson 


[34]).  In  addition,  since  the  matrices  (n)  and  G^ (n)  have  a  Toeplitz  structure 


the  vectors  (n)  am+i_j(n)  and  G^ (n)  am+i_j(n)  can  be  computed  by  using  fast  con 
volution  algorithms.  This  shows  that  the  number  of  operations  required  by  the  re 
cursions  (28)  and  (29)  is  0(mn  log  n) . 


i 


Increase  in  n 

To  compute  a(n+l,  m) ,  we  can  use  the  Toeplitz  structure  of  x.  To  do  so,  we 
consider  the  identity 


ao(n) 


x(n,m) 


Vn) 


\ 


\ 


Vi(n) 

\ 

\ 

\ 

am(n) 

Vi(n) 

a0(n) 

a(n,m) 

a(n,m-l) 

a(n,0) 

«(n) 


6(n)  0 

\ 

\ 

\ 

\ 

6(n) 


(30) 


where 


5  (n)  A 


0 

1 

| 

j 

I 

1 

0 

5  i 

j 

n  / 

n  +  1 


Then,  in  order  to  replace  the  row  by  row  scanning  of  the  set  of  D(n,m)  by  a  column 
by  column  scanning,  we  can  define  the  transformation 


P(n,m)  =  ®  In,j 


where  A  ®  B  denotes  the  Paley  product  of  two  matrices.  If  A  and  B  are  some  matrices 
of  size  n  x  n  and  ra  x  m,  this  product  is  defined  as 


A  consequence  of  the  structure  (27)  for  a(n,m)  is  that  the  residual  6^  does 
not  depend  on  m.  ,m 


A  0  B, 


A  ®  B  =  A  0  B. 

\  J 
\  A  0  B 


where  B^.  is  the  j  row  of  B.  By  multiplying  the  rows  and  columns  of  T(n,m)  by 
P(n,m),  and  by  observing  that  P(n,m)  =  PT(n,m)  =  P'^n.m),  one  can  rewrite  (30)  as 


T(n,m)  A(n,m)  » 


A(n,ra) 


Here,  we  have 


W-> 


Tj  (m) 


T(n,m)  = 


(32a) 


T9nj-1  0®) 


Wm) 


Ti(m) 


gil  gi0 

\  \ 


(32b) 


N  \ 

\  \ 

\  N 

gil  gi0 


? 


an(i,0) 


A0(n,m) 


A(n,m)  =  A^n.m) 


An(n,m) 


(33a) 


Ai (n,m) 


V1'1* 


an(i,m) 


a  (i,0) 
\ 


C33b) 


an ( i , 1 D  an(i,0) 


where  a  (i,j)  3  a  (i,j)  is  the  (i,j)th  coefficient  of  a  (z,w),  i.e. 
n  n,m  n  t  id 


n  m 


‘n,»  (Z’U)  "  I  I 


...  n-i  m-j 
an,m(l'j)  z  “ 


i»0  j-0 


(observe  that  the  structure  (27)  for  a(n,m)  implies  that  an  m(i»j)  does  not 
depend  on  m) .  Also,  the  m  +  1  *  ra  +  1  matrix  A(n,m)  is  given  by 

A(n,m)  *  diag  (6n>  .  (34) 

The  main  aspect  of  (31)  is  that  t(n,m)  and  A(n,m)  have  both  lower  triangular 
Toeplitz  (LT)  block  entries.  Since  the  algebra  of  LT  matrices  is  closed  (the 
product  of  two  LT  matrices  is  commutative  and  LT,  the  inverse  of  a  nonsingular 
LT  matrix  is  LT) ,  one  can  operate  on  the  blocks  and  as  if  they  were  scalars . 

This  means  that  to  solve  (31)  efficiently,  we  will  need  only  to  exploit  the  block 

A 

Toeplitz  structure  of  T(n,m).  This  will  be  done  by  using  a  set  of  recursions  de¬ 
rived  originally  by  Lanczos  [4],  and  introduced  in  the  context  of  realization  theory 


by  Kung  [35]  and  Kailath  [36]. 


These  recursions  are  based  on  the  observation  that  we  have 


where  R(n,m)  and  S(n,m)  are  LT  matrices  since  they  are  obtained  as  the  combina¬ 
tion  of  LT  matrices.  Consequently,  if  we  define 


M(n,m)  *  A'1  (n-1,  m)  A(n,m)  (36a) 

N(n,m)  =  A  *  (n,m)  R(n,m)  -  A  1  (n-1,  m)  R(n-1,  m)  ,  (36b) 


the  matrix 


A(n+1,  m) 


will  satisfy  the  equation  (31),  provided  that  we  replace  n  by  n  +  1.  In  this 
case,  the  new  residual  is  given  by 

A(n+1,  m)  *  S(n,m)  -  R(n,ra)  N(n,m)  -  S(n-1,  m)  M(n,m)  . 


(38) 


Since  LT  matrices  form  a  closed  algebra,  the  recursions  (37)  involve  only  the 
multiplication  ofm+ixm+ILT  matrices.  Therefore,  the  recursions  (37) 
require  only  0(nm  log  m)  operations  if  one  used  fast  Fourier  transform  techniques 
to  multiply  LT  matrices. 

Remark  The  LT  matrix  A(n+1,  m)  obtained  in  (38)  is  not  diagonal  in  general. 

A 

If  one  wants  A(n+1,  m)  to  be  diagonal  as  in  (34),  one  needs  only  to  factor 

A 

A(n+1,  m)  in  its  lower  triangular  part  with  unit  diagonal,  times  its  diagonal 

A 

part.  Then,  we  can  renormalize  A(n+1,  m)  accordingly. 

Another  useful  observation  is  that  the  residual  matrices  A(n,m)  =  diag  {6  } 

n 

have  to  be  assumed  nonsingular  for  the  previous  algorithm  to  be  valid.  When 
Sn  =  0  for  some  n,  a  generalized  set  of  recursions  (the  Berlekamp-Massey  recur¬ 
sions)  have  to  be  used.  These  recursions  were  introduced  in  the  1-D  scalar  case 
by  Berlekamp  [5]  and  Massey  [6] ,  and  then  generalized  to  the  1-D  matrix  case  by 
Dickinson,  Morf  and  Kailath  [37] .  However,  we  will  not  consider  this  case  here, 
and  it  will  be  assumed  that  6  j*  0  for  all  n. 


/ 
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III.  A  Modified  Approximation  Problem 

The  geometry  A(n,m;  2n+l,  m)  is  not  the  only  one  for  which  the  modified  Pade 
approximation  problem  (8)  has  a  recursive  solution.  In  this  section,  we  shall 
consider  the  geometry  A(n,m)  described  in  Figure  1.  As  mentioned  earlier,  the 
modified  approximation  problem  associated  with  A(n,m)  does  not  admit  an  equiva¬ 
lent  formulation  in  terms  of  rational  approximants .  However,  if  g^  m(z,w)  = 
b  (z,oj)/a  (z,io)  is  the  modified  Pade  approximant  associated  with  A(n,m),  we 
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shall  show  that  the  polynomials  a^  m(z,w)  and  b^  m(z,aO  obey  a  set  of  recursions 
similar  to  those  obtained  by  Jackson  [8]  for  2-D  orthogonal  polynomials  on  the 
hyper  real  line. 

/ 

The  first  step  is  to  decompose  A(n,m)  as  A(n,m)  =  -D(n,m)  UD+(n,m)  with 
D+(n,m)  =  { (i, j)  :  (i-1,  j-1)  e  D(n,m)}. 


Then,  by  considering  the  equation 


II.  .  .a  (z,cu)  g(z,0))  -  b  (z,w)  =  0  (39) 

A(n,m)  n,m  ’  '  '  n,nr  '  * 

we  can  verify  that  bn  m(z,w)  obeys  the  same  equation  as  in  Case  C  of  Section  II 
(it  is  given  by  (20)  or  (23)).  However,  an  m(z,ii))  satisfies  a  different  relation 
given  by 


V (n,m) 


a  (z,u>)  g(z,w) 
n,m 


=  <5  z 
n,m 


(n+1)  - (m+1) 


0) 


(40) 


provided  that  we  assume  that  the  coefficient  of  z~^n+^  u)  is  not  matched. 

Note  that,  as  in  Case  C  of  Section  II,  the  polynomial  a^  m(z,u>)  has  only 
(n+1) (m+1)  -  1  free  coefficients,  or  one  less  than  the  number  of  parameters  in 
D*(n,m).  If  a(n,m)  denotes  the  vector  obtained  by  scanning  the  coefficients  of 


an  m(z,w)  row  by  row,  (40)  can  be  tranformed  into 


T(n,m)  a(n,m)  =  e(n,m) 


where 


e(n,m)  * 


m  +  1  ® 


n  +  1 


Tj  Cn) 


T(n,m)  = 


(43a) 


T  ,  (n) 
m+1^ 


^+1) 


Tj  (n) 


(43b) 


s2n+lj 


The  matrix  T(n,m)  is  a  Toeplitz  block  Toeplitz  matrix,  and  it  corresponds  to  a 
simple  reordering  of  the  2-D  Hankel  matrix  H(n,m)  associated  with  g(z,w).  Indeed, 


if  one  denotes  by  the  k  x  k  matrix  given  by 


0  1 
/ 

/ 


/ 

1  0 


and  if  one  defines 


H(n,m)  =  T(n,m)  (Jffl+1  ®  Jfl+1), 


(44) 


/ 
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the  matrix  H(n,m)  has  a  Hankel  block  Hankel  structure.  It  differs  however  with 
the  Hankel  matrices  introduced  by  Fornasini  [13]  and  Kao  and  Chen  [15]  to  study 
the  2-D  realization  problem.  This  difference  will  be  explained  in  Section  IV, 
where  it  will  be  shown  that  the  Hankel  matrix  H(n,m)  does  not  play  as  central 
a  role  in  the  characterization  of  rational  functions  g(z,w)  =  b(z,u>)/a(z,w)  as 
it  did  in  the  1-D  case. 

The  Hankel  structure  of  H(n,m]  can  be  exploited  to  relate  the  partial  reali 
zation  problem  (41)  with  the  theory  of  orthogonal  polynomials  on  the  hyper  real 

2 

line.  To  do  so,  we  consider  a  nonnegative  weight  function  y(x,y)  defined  on  R 
and  we  denote  the  moments  of  y  by 


h.  .  = 


Then,  the  matrix 


/GO  v*  OO 

/  x1  y1  y(x, 

—00  J  -CO 


y)dxdy 


with 


H(n,m)  = 


Hj(n)  = 


J  -CO 

H0Cn) 

HjOO 

Hx(n) 

/ 

/ 

/ 

/ 

Hm^  ' 

hoi 

hu 

H  (n) 
m 


/ 


H2m^ 


h  . 


hij 


/ 


/ 


h  . 
nj 


h_  . 
2nj 


(45) 


(46a) 


(46b) 


can  be  used  to  characterize  the  inner  product  of  polynomials  of  degree  less  than 

n  and  m  in  x  and  y.  If  a(x,y)  and  b(x,y)  are  two  such  polynomials  and  if  we  de- 

*  * 

note  by  a  and  b  the  vectors  obtained  by  reversed  row  by  row  scanning  of  the 


coefficients  of  a(x,y)  and  b(x,y),  we  have 


<a,b>  = 


f  f 


It 

a(x,y)  b(x,y)  y(x,y)  dxdy  =  b  H(n,m)a  .  (47) 


Also,  if  we  consider 


n  m 


i=0  j  =0 


n-i  m-i 
a. .  x  y 


and  if  a  and  a  denote  the  vectors  of  coefficients  of  a(x,y)  obtained  by  direct 
and  reversed  row  by  row  scanning,  one  has 


a  =  (J  .  ®  J  ,)a 
m+1  n+1 


where 


a  =  la 


Consequently,  if  a^  mCx,y)  is  a  monic  polynomial  such  that  an  mtx>>0  ±  xXyJ  for 
(0,0)  <  (i,j)  <  (n,m)  and  if  one  denotes 


<a  ,  xnyI11>  »  5 

n,m  n,m 


the  vector  a  (n,m)  of  coefficients  of  a^  satisfies  the  equation 


H(n,m)  a  (n,m)  =  e(n,m) 


where  e(n,m)  is  defined  as  in  (42).  By  taking  (44)  and  (48)  into  account,  this  | 

shows  that  the  coefficients  of  an  m(x»y)  satisfy  exactly  the  same  equation  as  the  j 

coefficients  of  a  (z,oj)  in  the  partial  realization  problem.  This  observation  sug-| 
n,m  ] 

gests  that  the  recursions  derived  by  Jackson  [8]  for  the  orthogonal  polynomials  J 


I 
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an  m(x,y)  can  be  adapted  to  the  approximation  problem  considered  here.  The  only 

difference  between  these  two  problems  is  that  in  the  study  of  2-D  orthogonal 
2 

polynomials  on  iR  ,  the  Hankel  matrix  H(n,m)  of  the  moments  of  u(x,y)  is  nonnegative 
definite,  while  this  is  not  usually  the  case  when  H(n,m)  is  constructed  from  the 
Markov  parameters  of  g(z,a>).  However,  as  we  shall  see,  this  difference  does  not 
play  a  role  in  the  recursions  satisfied  by  an  m(z,u)). 


Auxiliary  Solutions 

One  of  the  main  features  of  Jackson's  recursions  and  of  the  recursions  that 
we  shall  present  below  is  that  they  require  the  introduction  of  several  auxiliary 
solutions.  At  stage  (n,m),  we  will  introduce  n  +  m  +  1  modified  Pade  approximants 
which  will  be  divided  into  n  +  1  horizontal  approximants  and  m  +  1  vertical  approxi¬ 
mants.  The  horizontal  approximants  will  be  denoted 


g(z,«*»)  =  k*  (z,a))/h*  (z,u>) 

n,m  n,m 


with  0  _<  i  £  n.  They  are  such  that 

n 

.  i  ,  ,  i  m  \  ,  i  ..  ...  n-k  m 

h  (z,u>)  =  z  a)  +  /  h  (k,0)z  w 


n,m 


k=n-i+l 


n,m 


(50) 


n  m 


2  2 


-k  m-Z 

OJ 


k=0  £=1 


(51) 


a  i  i 

is  an  asymetric  South-west  (Sw)  causal  filter  and  such  that  g  =  k  /h  is  a 
7  v  J  6  n,m  n,m 

modified  Pade  approximant  of  g  over  the  domain  A  described  in  Figure  3a.  This 
means  that  k*  m(z,uj)  and  h*  m(z,uj)  satisfy,  respectively,  the  equations 


V(n,m)  0Z'“)  *tz’“>  =  S-  -(z,z 


-(n+1)  -(m+1) 


n,m 


u> 


(52a) 

(52b) 
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where  61  (z)  is  a  polynomial  in  z  only  such  that 

n,ra 


deg  (z)  =  n-i. 
n,m 


C53) 


Similarly,  the  vertical  approximants  will  be  denoted 


g(z,w)  =  (z (z,u>) 
n.  >  in  n,m 


(54) 


where 


m 


(z  ,u>)  =  z 

n,m 


v  *  2 
£=m- j+1 


m-Z 


♦2  i 

k=l  £=0 


n-k  m-Z 


(55) 


is  a  West-south  (Ws)  filter  with  0  <  j  <  m  and  where  g  =  uJ  is  a  modified 

—  —  n,m  n,m 

Pade  approximant  of  g  for  the  geometry  A  described  in  Figure  3b.  Therefore,  the 

polynomials  U-*  (z,<d)  and  vJ  (z ,w)  satisfy  the  equations 

n,m  n,m 


IIn+  r  (Z.OJ)  g(z,u)  =  (oo)  z 

D+(n,m)  n,nr  ’  3  '  'n,nr  ‘ 


(56a) 

(56b) 


where  Y"1  (w)  is  a  polynomial  in  w  only  such  that 

n  f  ni 


degY^,m(u))  =  m-j. 


(57) 


The  total  number  of  auxiliary  approximants  is  only  n  +  m  +  1  if  one  observes  that 

g(z,u>)  =  k”  fz,u))/h"  (z,cj)  *  u™  (z,u)/v"  (z,w) 

n,m  n,m  n,m  n,m 


=  b  (z,aj)/a  (z,co) 
n,nr  ’  n,m  ’ 


(58) 
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is  the  quarter  plane  Pade  approximant  that  we  want  to  compute. 

The  equations  (52)  and  (56)  can  be  written  in  matrix  form  by  scanning 

the  coefficients  of  kj  ,  h*  ,  u*'  and  row  by  row  and  by  denoting  by 

n,m  n,m  n,m  n,m  7  &  / 

K(n,m)  =  (k°(n,m),  ....  kn(n,m)) 


L(n,m)  *  (h°(n,m),  ....  hn(n,m)) 

U(n,m)  =  (u°(n,m),  ....  um(n,m)) 
V(n,m)  =  (v°(n,m),  ....  vm(n,m)) 


the  block  matrices  obtained  by  grouping  the  vectors  of  coefficients  of  k1  ,  h1 

n,m’  n,m 

and  v?  ,  .  Then,  the  numerators  satisfy  the  relation 

n,m  n,m  ' 


(K(n,m),  U(n,m))  =  G(n,m)  (L(n,m),  V(n,m)) 


(59) 


where  G(n,m)  is  given  by  (15),  so  that  the  main  problem  is  to  compute  l(n,m)  and 
V(n,m).  To  do  so,  we  shall  use  the  identity 


T(n,m)  (L(n,m) ,  V(n,m))  =  (D(n,m),  C(n,ra)) 


(60) 


where  T(n,m)  is  given  by  (43)  and  D  and  C  are  matrices  of  respective  size 
(n+1) (m+1)  x  (n+1)  and  (n+l)(ra+l)  x  (m+1)  which  are  given  by 

I  ■' 


D(n,ra)  = 


0 

i 

i 

I 

I 

0 


m+1  ®  A(n,m) 


(61a) 


and 


/ 
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n  +  1  . 


(61b) 


Here,  A(n,m)  and  r(n,m)  are  the  n  +  1  x  n  +  1  and  m  +  1  x  m  +  1  matrices  whose 

columns  A1 (n,m)  and  rJ(n,m)  are  given,  respectively,  by  the  coefficients  of  the 

residual  polynomials  61  (z)  and  (to).  An  important  feature  of  A(n,m)  and 

n  f  in  n  jin 

r(n,m)  is  that  these  matrices  are  lower  triangular.  This  property  is  a  consequence 
of  (53)  and  (57). 


Increase  in  m 

The  introduction  of  the  auxiliary  solutions  L(n,m)  and  V(n,m)  is  motivated 
mainly  by  the  fact  that  these  solutions  can  be  computed  recursively.  To  do  so, 
we  shall  use  an  algorithm  based  on  the  block  form  of  the  recursions  discussed  in 
Lanczos  [4].  The  main  difference  with  Lanczos'  recursions  is  that  when  m  is  in¬ 
creased  to  m  +  1,  one  has  not  only  to  compute  L(n,m+1),  but  also  V(n,m+1)  (note 
that  this  requires  the  construction  of  one  additional  auxiliary  solution  since  at 
stage  (n,m+l),  there  are  m  +  2  vertical  approximants) .  Thus,  we  have 


autation  of  L(n,m+lj 


One  uses  the  identity 


T(n,m+1) 
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0 


(62) 


which  is  similar  to  (35)  with  one  important  modification:  instead  of  being  lower 
triangular  Toeplitz  (as  for  T(n,m))  the  blocks  of  T(n,m)  are  fully  Toeplitz.  This 
difference  is  significant  since  the  algebra  of  Toeplitz  matrices  (unlike  the  alge¬ 
bra  of  lower  triangular  Toeplitz  matrices)  is  not  closed.  In  fact,  its  closure 
is  formed  by  the  sums  of  products  of  lower  times  upper  Toeplitz  matrices  (see, 
e.g.,  [38],  [39]).  This  means  that  the  block  entries  of  L(n,m)  as  well  as  A(n,m), 
R(n,m),  S(n,m)  are  not  Toeplitz  in  general.  These  matrices  are  usually  arbitrary. 
Thus,  if  as  in  (36),  we  define 

M(n,m)  =  A  *(n,m-l)  A(n,m)  (63a) 

and 

N(n,m)  =  A  ^n.m)  R(n,m)  -  A"*(n,m-1)  R(n,m-1)  (63b) 

the  recursions 


L(n,ra+1) 


1 

/  °  1 

/  ° 

. 

L(n,m) 

- 

N(n,m)  - 

L(n,m) 

L(n,m-1) 

o 

M(n,m) 


(64) 


/ 
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require  0(mn)3  operations  instead  of  0(mn  log  n)  for  the  corresponding  recursions 
in  the  geometry  C  of  Section  II. 

Computation  of  V(n,m+1) 

The  first  step  is  to  observe  that 


1  °  1 

T(n,m+1) 

1 

L(n,m+1)  = 

C(n,m) 

0 

V(n,m) 

I 

P(n,m) 

A(n,m+1) 

(65) 


Therefore,  if  we  introduce 

Q(n,m)  =  A  *(n,m+l)  P(n,ra)  (66a) 

and 


/ _ o 


V(n,m+1) 


\  V(n,m) 


-  L(n,m+1)  Q(n,m)  , 


(66b) 


the  matrix  V(n,m+1)  satisfies  the  equation 


T(n,m+1)  V(n,m+1) 


C(n,m) 


(67) 


0 
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This  implies  that  the  first  m  +  1  columns  of  V(n,m+1)  are  given  by  V(n,m+1). 

To  obtain  the  new  auxiliary  solution  (i.e.  the  last  column  of  V) ,  one  needs  only 
to  note  from  (58)  that  the  last  column  of  V(n,m+1)  is  the  same  as  the  last  column 
of  L(n,m*l).  But  this  column  has  just  been  computed,  so  that 

V(n,m+1)  *  (V(n,m+1),  hn(n,m+l‘J) .  (68) 

The  number  of  operations  required  by  (68)  is  also  O(mn^). 


Increase  in  n 

Due  to  the  symmetry  of  the  domain  of  approximation  A(n,m),  to  increase  n 
one  can  use  the  same  recursions  as  for  m,  provided  that  the  set  D(n,m)  is  scanned 
column  by  column  instead  of  row  by  row,  and  that  the  roles  of  L(n,m)  and  V(n,m) 
are  exchanged,  as  well  as  those  of  n  and  m. 

Remark  1  There  is  a  significant  difference  between  the  previous  recursions  for 

the  Pade  approximants  h^  ^  and  v^  and  those  derived  by  Jackson  [8]  for  2-D 

2 

orthogonal  polynomials  on  IR  .  This  difference  arises  from  the  fact  that  there  are 
several  complete  orderings  of  the  monomials  x1  'p  which  are  compatible  with  the 
partial  ordering 

x1  p  <  xk  yl  iff  i  <  k,  j  <  l. 


using  the  polynomials  of  degree  n  and  n  -  1  to  compute  those  of  degree  n. 

By  comparison,  the  method  that  we  have  used  to  compute  h*  m  (resp.  m)  is 

based  on  a  simple  row  by  row  (resp.  column  by  column)  lexicographic  order  of  the 

monomials  z1  uP .  In  this  framework,  instead  of  performing  a  complete  row  by  row 
2 

scan  of  Hi  (the  rows  would  be  infinite),  we  truncate  the  rows  to  length  n,  i.e., 
we  consider  the  monomials  z1  such  that  i  <  n,  and  we  have 


i  j  k  t 

Z  tii  <  Z  OJ 


if  either  i  <  £  or  j  =  £  and  i  <  k.  Then,  the  recursions  for  h1  involve  either 
J  J  n,m 

an  order  increase  (m  -*•  m+1)  or  a  change  of  truncation  (n  -»■  n+1). 


Remark  2  Since  a(n,m)  *  hn(n,m)  =  vm(n,m),  the  previous  recursions  enable  us 
to  compute  a(n,m) .  To  compute  b(n,m),  one  needs  only  to  observe  that  b(n,m)  = 
kn(n,m)  *  ura(n,m)  and  that 

(K(n,m),  U(n,m))  =  G(n,m)  (L(n,m),  V(n,m)) 


so  that  the  matrices  K(n,m)  and  U(n,ra)  obey  exactly  the  same  recursions  as  L(n,m) 
and  V(n,m).  The  same  observation  holds  for  the  computation  of  the  numerator 
b(n,m)  in  the  geometry  C  of  Section  II. 


IV.  Complete  Realization  of  2-D  Transfer  Functions 


The  previous  partial  realization  schemes  can  be  related  to  the  results 
obtained  by  Fomasini  [13],  Clerget  [14]  and  Kao  and  Chen  [IS]  for  the  study 
of  the  complete  realization  problem.  In  this  case,  we  want  to  characterize  the 
2-D  quarter-plane  transfer  functions  g(z,u>)  which  are  rational,  i.e.,  such  that 


g(z,w)  *  b(z,aO/a(z,oO 


and  such  that  a(z,w)  has  degree  (n,m). 

Then,  if  one  considers  the  geometry  B  of  Figure  2,  for  all  (p,q)  >,  (n,m), 
the  rational  function  g(z,i»0  satisfies 

V,m;p,q)  8(2>“°  '  l’U>“)  "  °  (69) 

where  A(n,m;p,q)  is  defined  as  in  (10).  This  means  that  by  decomposing  the  domain 
A(n,m;p,q)  as  A(n,m;p,q)  »  -D(n,m)  uA(n,m;p,q),  where  A(n,m;p,q)  is  the  complement 
of  -D(n,m)  in  A(n,m;p,q),  one  has 


V,n,;p,q)  3(z-“)  =  0 


This  identity  can  be  used  to  characterize  the  transfer  functions  g(z,co)  which  admit 
a  rational  realization  of  degree  (n,m).  To  do  so,  we  rewrite  (70)  in  matrix  form 
by  introducing  the  matrices 


T(n,m;p,q) 


Tm+i (n;p)  Tm(n;p)  •  '  •  Ti (n^P) 
Tm+2(n;p:)  Tm+l(n;p)-  *  •  T2(n;p) 


(71a) 


T.(n;p)  T_  ^njp).  .  .  Tqm(n;p) 


/ 

•*»*-*:  *31 W*  y.-  , 
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T(n,m;p)  = 


T0(n;p) 

TjCnjp)  TQ(n;p) 


|  Tm(n;P) 


T1  Cn;p)  TQCn;pl 


0(n,m;q)  = 


where 


Tj (n;p) 


W") 

em  (n)  .... 

^(n) 

0m+2™ 

9ra+l(n)  .... 

02(n) 

eq(n) 

0q.l(n)  .... 

e  "  o 

q-m^ 

8n+lj 

gnj  *  *  *  * 

glj 

gn+2j 

gn+lj  ■  •  •  • 

g2j 

®Pj 


*P-lj  ’ 


P-nj 


and 


ej (n;p) 


s0j 


nJ 


g0j 


g0j 


(71b) 


(71c) 


Then,  if  we  define 


/ 

* 
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I (n,m;p,q)  * 


/  x (n,m;p) 
T(n,ra;p,q) 
9(n,m;q) 


(72) 


and  if  a  denotes  the  vector  obtained  by  scanning  the  coefficients  of  a(z,uj)  row 
by  row,  (70)  can  be  written  as 

Z(n,m;p,q)a  =  0  .  (73) 

This  identity  leads  to  the  following  realization  criterion  for  g(z,w). 

Theorem  (cf.  Kao  and  Chen  [15]) 

If  E(n,m)^  E(n,m;  2n+l,  2m+l),  the  transfer  function  g(z,w)  has  a  rational 
realization  of  order  (n,m)  if  and  only  if 

rkE(n,m)  =  rkZ(n,m;p,q)  =  (n+l)(m+l)  -  1  (74) 

for  all  (p,q)  (2n+l,  2m+l),  and  if  the  first  nm  +  n  +  m  columns  of  Z(n,m) 

are  independent. 

The  main  aspect  of  this  result  is  that  it  depends  on  the  3(n+l) (m+1) x(n+l) (m+1) 
matrix  E(n,m)  which  does  not  have  either  a  Toeplitz  or  a  Hankel  structure  (note, 
however,  that  in  Kao  and  Chen  [15],  this  matrix  is  said  to  be  the  Hankel  matrix 
of  g(z,w)).  By  identifying 


x(n,m)  =  r(n,m;  2n+l) 

T(n,m)  =  T(n,m;  2n+l,  2n+l) 


(75) 


with  the  matrices  defined  in  (24)  and  (43),  it  is  clear  that  the  matrix  T(n,m) 
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(which  is  obtained  by  reordering  the  columns  of  the  2-D  Hankel  matrix  H(n,m )) 
is  only  a  submatrix  of  Z(n,m).  This  shows  that,  unlike  in  the  1-D  case, 
H(n,m)  cannot  be  used  to  characterize  the  rationality  of  g. 


/ 

* 
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V.  Stochastic  Modeling 

In  this  section,  the  2-D  realization  problem  will  be  considered  from  a 

stochastic  point  of  view.  It  will  be  assumed  that  we  are  given  a  2-D  zero-mean 

2 

stochastic  process  y(i,j),  (i,j)  e  2  ,  which  is  space- invariant,  so  that  its 
covariance  is  such  that 

Tjj  ■  E[y(k+i,  A+j)y(k,4)]  (76) 

2 

for  all  (k,£)  e  2  .  Then,  we  shall  seek  to  find  a  model  for  y(i,j)  which  is 
both  autoregressive  and  causal,  so  that 

y(ij)  +  J  ak^y(i-k,  j-A)  =  u(i,j)  (77) 

I- (0,0) 

where  u(i,j)  is  a  2-D  white  noise,  i.e., 

E  [u(i,  j)u(k,A)  ]  =  u6i_k6^_Jl  , 

2 

and  where  I  is  a  finite  and  causal  subset  of  2  (by  causal,  we  mean  that  I  be- 

2 

longs  to  an  asymmetric  half  plane  of  2  ) .  The  main  property  of  such  models  is 
that  if 

r(z,u0  =  ^  rijZ~ia)"^ 

is  the  spectrum  of  y(*,*)>  and  if 

a(z,tu)  =  1  +  ^  a^z-1^"3 

I-(0,0) 


is  the  polynomial  associated  with  the  filter  (77),  then 


/ 
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r(z,w)  = 


u 


a(z,  (lijafz"1,^"1) 


(78) 


is  a  spectral  factorization  of  r(z,w).  The  2-D  spectral  factorization  problem 
has  already  been  the  object  of  a  large  amount  of  attention  (cf.  Helson  and 
Lowdenslager  [28],  Pistor  [29],  Ekstrom  and  Woods  [16],  Murray  [30],  Delsarte, 
Genin  and  Kamp  [27]),  and  here  this  problem  will  be  considered  only  indirectly. 
However,  one  needs  to  recall  the  conditions  of  existence  of  factorizations 
such  as  (78).  If  I  is  chosen  to  be  an  asymmetric  half-plance,  i.e., 

(i)  (0,0)  el 

(ii)  (i,j)  £  I  if  and  only  if  (-i,-j)  t-  I,  unless  (i,j)  =  (0,0) 

(iii)  if  (i, j)  and  (i',j')  e  I,  then  (i+j ,  i'+j')  e  I, 
and  if 

/log  r(e  ,e1<P)d6d<j>  >-»  ,  (79) 

4tt  0  J0 

it  is  shown  in  [28]  that  there  always  exists  a  factorization  (78)  such  that 
a(*,*)  is  stable.  However,  the  support  I  of  a(*,*)  is  in  general  infinite. 

If  I  =  IN  is  the  positive  quarter  plane,  it  is  shown  in  Murray  [30]  that  the 
factorization  does  not  exist,  unless  the  Fourier  coefficients 

2 it  r2ir 


kA 


=  f  [  log  r(e:*'®,  e*^)  exp.  i  (k6  +  £<J>)d0d<{> 

Jr\  Jf\ 


4tt  0  '0 


are  zero  in  the  second  and  fourth  quadrant  of  2  ,  i.e.,  s^  =  0  if  k£  <  0. 
In  this  section,  we  will  be  concerned  with  the  problem  of  finding  some 


approximants  for  the  spectral  factor  a(*,*)  .  This  approximation  problem  will 
be  formulated  as  a  causal  linear  least-squares  estimation  problem  for  the  process 
y(*»*).  At  stage  (n,m),  we  shall  consider  the  set  of  asymmetric  half-plane 


7 


m 
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predictors  of  order  (n,m)  for  y  (*,*)•  There  are  n+m+1  such  predictors  and, 

as  in  Section  III,  they  can  be  divided  into  n+1  horizontal  predictors  h1  0<i<n, 

n,m, - 

and  into  m+1  vertical  predictors  v^  m»0<p<m.  The  horizontal  predictors  are 
such  that 

n 

k=i+l 


n  m 


•  2  X  hi,mCk’*)yCt-k’  s-£) 


(80) 


k=0  4=1 


is  the  linear  least-squares  estimate  of  y(t-i,  s)  given  the  observation  of  y(*,«) 
over  the  set  H*  m  described  in  Figure  4. a  (where  0<i<n) .  These  predictors  have 
the  property  that  the  error 


e"-i-  sK,m>  ■  s>  -  i’!*-1- 

is  orthogonal  to  y(*,»)  over  the  domain  H1  ,  i.e. 

n,m 

E[e(t-i,  s|HjJ  m)y(p,q)]  =  0 

for  all  (p,q)  eH^.  By  introducing  the  vector 

yn (t , s ;  n)  \ 


Y(t,s;  n,m)  =  |  y..(t,s;  n) 
yra(t,S;  n)/ 


(81) 


(82) 


(83a) 


where 


y(t,  s-j) 


(83b) 


yjCt.s;  n) 


I 

,  y(t-i,  s-j) 

\  y(t-n,  s-j) / 


is  the  jt^  row  of  data  in  ffl,  and  by  noting  that 

e(t-i,  s|H*  )  =  YT(t,s;  n.nOhV.m)  (84) 

a,ra 

where  h1 (n,m)  is  the  vector  obtained  by  scanning  the  rows  of  the  filter  h*  m(*,*), 
one  can  rewrite  (82)  as 

R(n,m)h* (n,m)  *  e^(n,m)  (85) 

T 

where  R(n,m)  =  E[Y(t,s;  n,ra)Y  (t,s;  n,m)]  is  the  covariance  of  Y(t,s;  n,m), 
and 

1 
0 

0 

Here,  61(n,m)  is  the  n+l-vector  obtained  by  correlating  e(t-i,  s|H*  m)  with  the 
row  yQ(t,s;  n) .  This  means  that  the  last  n-i  entries  of  61(n,m)  are  zero,  but 
fi'in.ffl)  itself  is  always  nonzero  if  the  prediction  problem  is  nonsingular,  i.e., 
if  the  mean-square  prediction  error  E[e^(t-i,  s|H^  m)]  ^0.  The  equation  (84) 
presents  the  feature  that  R(n,m)  is  Toeplitz  block  Toeplitz.  Indeed,  we  have 


ra+1  0 


(86) 
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R(n,m)  = 


R0(n) 

R^n).  Rm(n) 

• 

R.l(n) 

R0(n)  *  • 

Rj  (n) 

. 

R-m(n) 

•  R_x (n)  *  R0(n) 

(87a) 


with 


RjCn)  = 


roj 


-lj  . 


-n] 


roj 


-lj 


nj 


*r. 


•  r 

r0j 


(87b) 


Our  main  objective  here  will  be  to  take  advantage  of  this  structure  to  obtain 

some  recursions  for  the  horizontal  and  vertical  predictors  h1  and  v^ 

r  n,m  n,m 

The  vertical  predictors  are  defined  symmetrically.  Thus, 

m 

=  -  2  vi,m(0’*)yCt’  S-A) 

i-i +i 


m 


-2  2  vj>m(k,*)y(t-k,  »-« 


k=l  z=o 


(88) 


is  the  linear  least-squares  estimate  of  y(t,  s-j)  given  the  observation  of 
y(*j*)  over  the  set  m  described  in  Figure  4.b  (where  0<j£m) .  In  this  case, 

e(t,  s-jlvjn,m)  =  y(t>  s-V  '  s-jl^n.m5  (89) 

is  orthogonal  to  y(*,*)  over  the  domain  so  that  by  denoting  by  v^(n,m) 


firrmunm- 


/ 
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the  vector  obtained  by  scanning  the  coefficients  of  vJ  (•,*)  row  by  row  we 

n,m  1  ' 

obtain 


R(n,m)vJ  (n,m)  »  n3(n,m) 


(90) 


The  vector  (n,m)  is  given  by 


n^fn.m)  *YJ(n,m)©  n+1 


(91) 


where  yJ (n,m)  is  the  m+l-vector  obtained  by  correlating  e(t,  s- j | )  with  the 

n.y 

column  vector 


(y(t,  s)  ' 

X(t,  s-j)  . 
y(t.  s-m) 

By  noting  that 

h°(n,m)  »  v°(n,m)  =  a(n,m)  (92) 

where  a(n,m)  is  the  quarter-plane  predictor  of  order  (n,m)  associated  with  y (*,*)» 
we  see  that  the  number  of  predictors  introduced  at  stage  (n,m)  is  only  n+m+1. 

Then,  if  one  denotes 


H(n,m)  *  (h®(n,m) . . .hn(n,ra)) 
V(n,m)  *  (v°(n,m) . . .vm(n,m))  , 


(t-n,s-m) 


Figure  4. a 


Figure  4.b 


The  Prediction  Geometry  for  the  Horizontal  Predictor  h1 

n, 


(t,s) 

(t,s-(j+0) 


The  Prediction  Geometry  for  the  Vertical  Predictor 


/ 
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the  equations  (85)  and  (90)  can  be  regrouped  as 


R(n,m)(H(n,m),  V(n,m))  =  (e(n,m),n(n,m))  (93) 

where  e(n,m)  and  q(n,m)  are  the  matrices  whose  columns  are  e^n.m)  and  nJ(n,m). 
To  solve  (93) ,  we  shall  now  introduce  a  2-D  generalization  of  the  recursions 
introduced  by  Levinson  (19],  and  by  Krein  [40]  in  the  continuous  case,  to  solve 
1-D  Toeplitz  equations. 

The  recursions  that  we  consider  were  first  presented  in  [20],  [21]  and 
their  main  feature  is  that  n  and  m  can  be  increased  separately.  By  comparison, 
we  note  that  the  recursions  introduced  by  Justice  [22]  were  requiring  that 
either  n  or  m  be  fixed  a  priori  (in  fact  Justice's  recursions  were  more  like 
those  introduced  by  Whittle  [41]  and  Wiggins  and  Robinson  [42]  to  generalize 
Levinson's  recursions  to  the  1-D  matrix  case).  These  recursions  differ  also 
from  those  considered  by  Marzetta  [24]  and  Delsarte,  Genin  and  Kamp  [18]  for 
asymmetric  half-plane  Toeplitz  systems  in  the  sense  that  the  recursions  derived 
by  these  authors  were  corresponding  to  a  different  geometry  involving  infinite 
vertical  (or  horizontal)  strips. 

The  first  step  in  the  derivation  of  the  2-D  Levinson  recursions  is  to  intro 
duce  the  filters 


which  are  obtained  by  reversing  the  direction  of  propagation  of  the  horizontal 

and  vertical  predictors  h1  and  vJ  .  These  filters  can  easily  be  seen  to  pro- 

n,m  n,m 


vide  the  linear  least-squares  North-east  (Ne)  and  East-north  (En)  asymmetric 
half-plane  predictors  of  order  (n,m)  associated  with  the  process  y  (*,•)■  Then, 

*i  *•; 

the  vectors  of  coefficients  of  h  and  v  J  are  given  by 

n,m  n,m 

h^Cn.m)  ■  (Jm+1®  Jn+1)h1(n,«)  - 

v*^  (n,m)  *  (Jn+1  ®  (n,m)  (95) 

♦  * 

and  if  one  denotes  by  H  (n,m)  and  V  (n,m)  the  matrices  obtained  by  regrouping 
these  vectors,  one  finds  that 

R(n,m)(H  (n,m),V  (n,m))  =  (e  (n,m),  n*(n,m))  (96) 


where 


* 

e  (n,m) 


0  |®A  (n,m),  A  (n,m)  =  Jn+^A(n,m) 


(97) 


n  (n,m)  *  r  (n,m)  ®  I  ^  |  *  r  (n»m)  a  Jm+irCn»m)> 


the  matrices  A(n,m)  and  r(n,m)  being  the  matrices  whose  columns  are  S^n.m)  and 
y^(n,m).  Now,  the  2-D  Levinson  recursions  can  be  described  as  follows: 


Increase  in  m 

To  increase  m,  one  has  not  only  to  compute  H(n,m+1)  (this  will  be  done  by 


using  the  block  1-D  Levinson  recursions),  but  also  V(n,  m+1),  a  step  that  in 
volves  the  introduction  of  one  auxiliary  solution. 


-45- 


Computation  of  H(n,  m»l) 
Consider  the  identity 


(98) 


where  the  residuals  are 

a(n,m)  =  (R_(m+1)(n)  R_1(n))HCn,m) 


0  | 

1  A(n,m) 

a*(n,m) 

R(n,  m+1) 

H(n,m) 

H  (n,m)  1 

s 

0 

0 

0 

1  a(n,m) 

A* (n,m) 

a  (n,m)  ^  Jn+1<*(n>m) 


Then,  if  we  define 


H(n,m+1)  = 


H(n,m) 


with 


H  (n,m) 


p(n,m) 


(99) 


*-l 

p(n,m)  »  A  (n,jn)oi(n,m)  ,  (100) 

the  matrix  H(n,m+1)  satisfies  (93),  where  we  have 

* 

A(n,m+1)  =  A(n,m)  -  a  (n,m)p(n,m)  .  (101) 

However,  A(n,m+1)  is  not  upper  triangular  in  general  as  would  be  required  by 


f 
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the  geometry  of  horizontal  prediction  chosen  here.  To  transform  it  into  this 
form,  one  needs  only  to  factor  A(n,m+1)  in  its  upper  times  lower  part  and 
renormalize  H(n,m+1)  accordingly. 


Computation  of  V(n,m+1) 

One  needs  to  note  that 


R(n,m+1) 


V(n,m) 


H(n,m+1) 


8  (n,m) 

A (n,m+l)  \ 

n  (n,m) 

• 

(102) 


where 


3 (n,m)  4  (Rx (n)  ...  Rm+1 (n))V(n,m) 
Thus,  if  we  introduce 

a(n,m)  =•  A-1  (n,m+l)S (n,m)  , 

the  matrix 


/  \ 

/ 

1  ° 

V(n,m+1)  »  1 

- 

H(n,m+1) 

\  V(n,m)  j 

l  / 

(103) 


(104) 


obeys  the  same  equation  as  the  last  m+1  columns  of  V(n,m+1).  To  obtain  the 
first  column,  one  needs  only  to  note  that  the  first  columns  of  V  and  H  are  iden¬ 
tical,  so  that 

V(n,m«-1)  *  (h°(n,m*l),  9(n,m+l)) 


(105) 
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In  general,  as  noted  earlier,  even  though  the  blocks  of  R(n,ra)  have  a  Toeplitz 
structure,  this  is  not  the  case  for  H(n,m)  and  V(n,m),  so  that  the  number  of 
operations  required  by  the  recursions  (99)  and  (104)  is  O(mn^).  Also,  a  condi¬ 
tion  for  the  validity  of  the  previous  recursions  is  that  A(n,m)  be  nonsingular. 
But,  since  A(n,m)  is  upper  triangular,  and  since  its  diagonal  terms  are  equal 
to  E[e^(t-i,  s|H^  m)J,  0<i<n,  we  see  that  A(n,m)  will  be  invertible  if  the 
estimation  problem  is  nonsingular. 

Increase  in  n 

Similar  recursions  can  be  obtained  by  reordering  R(n,m)  in  blocks  of  size 
m+1  x  m+1,  and  by  exchanging  the  roles  of  n  and  m  and  those  of  H  and  V. 

Remark .  The  previous  recursions  can  be  related  to  those  obtained  by  Genin  and 

Kamp  [23]  for  2-D  orthogonal  polynomials  on  the  unit  hypercircle.  This  relation 

is  based  on  the  isomorphism  [43]  existing  between  the  Hilbert  space  of  random 

2 

variables  y(i,j),  (i,j)  e  iN  with  scalar  product 

<y(i,j),  y».i>  -  rRH 

2 

and  the  space  of  functions  which  are  square  integrable  over  [0,2tt]  with  respect 

1  A  4  x 

to  the  positive  weight  function  r(e  ,e  *).  From  this  point  of  view,  the  predic¬ 
tion  problem  in  the  plane  and  the  one  of  orthogonalizing  the  monomials  z1^  with 
respect  to  r(»,»)  on  the  hypercircle  |z|  =  |u>|  =  1  are  exactly  the  same.  However, 
as  in  the  deterministic  realization  problem,  there  is  a  difference  between  our 
recursions  and  those  of  Genin  and  Kamp  which  is  due  to  the  total  order  of  the 
monomials  z1^  that  we  have  chosen.  As  in  the  deterministic  case,  we  have  cho¬ 
sen  a  lexicographic  order  with  truncation,  while  Genin  and  Kamp's  order  was 
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based  on  the  grade  of  the  monomials  z1^,  i.e., 
grade  z1^  =  max(i,j)  , 

so  that  one  had  z'V  <  zk<^  if  either  grade  z^oi1  <  grade  zko/  or  grade 
zV  *  grade  zki/  and  j  <  l  or  k  <  i. 

The  previous  recursions  enable  us  to  compute  the  planar  predictors  h*  and 

vJ  recursively.  However,  as  was  mentioned  at  the  beginning  of  this  section, 
n,m 

the  motivation  for  computing  these  predictors  is  to  approximate  the  spectral 

factors  of  r(z,<D).  We  now  justify  this  claim  by  considering  the  stability  and 

convergence  properties  of  h1  (z , toj  and  v*1  (z,oj).  The  first  observation  is 

n,m  n  j  m 

that,  unlike  in  the  1-D  case,  the  filters  h1  (z,w)  and  v ^  (z,w)  are  not  always 

n,m  n,m  - * — 

stable,  as  noted  in  [23],  unless  r(z,u>)  is  separable,  i.e.. 


r(z,w)  =  ^(zj^Cw). 

However,  it  was  shown  by  Helson  and  Lowdenslager  [28],  and  by  Delsarte,  Genin 
and  Kamp  [18],  that 


Theorem:  Convergence  of  h1 
_ _ _ n,m 

If  the  function  rCe1®^1^)  is  strictly  positive  and  summable  over  [0,2ir]2,  and 
if  agw(z,ui)  denotes  the  spectral  factor  (78)  associated  with  the  South-west  Half¬ 
plane  I  »  {(i,j):  j  <  0,  or  j  »  0  and  i  <  0}  ,  when  n-k,m  and  k  «,  we  have 


hk 

n,m 

2 

over  [0,2ir]  . 
Since  aSw 


(ei9,e^)  *  aSw(ei0,e^) 


(106) 


(z,w)  is  stable,  this  shows  that  when  n,m  and  i  tend  to  infinity. 


/ 


m(z,od)  is  stable.  In  fact,  it  was  shown  by  Chang  and  Aggarwal  [17]  that,  for 

a  fixed  value  of  m,  there  exists  an  integer  such  that  h1  (z,co)  is  stable  for  i > k„ 

0  n,m  —  0 

and  n-i  kQ.  By  symmetry,  when  n,m-j  -*■  80 ,  one  gets  m(ei8,ei(^)  -*■  a^fe  .e1^) 

where  aWg  is  the  spectral  factor  corresponding  to  the  West-south  half-plane. 

However,  these  results  do  not  settle  completely  the  convergence  problem 
for  the  filters  h1  or  v^  .  For  example,  the  existence  of  limits  for  h1  fz ,oj) 
when  n,m  -*■  00  with  i  constant  is  not  clear.  In  the  special  case  when  i  =  0, 
h8  m(z,u)  =  a^  is  the  quarter-plane  predictor  of  order  (n,m)  associated 

with  y(* ,  •)  and  only  partial  results  on  the  convergence  of  an  m(z,u))  are 
available  (see  [2S]-[27]).  In  fact,  this  question  is  related  to  Shanks'  con¬ 
jecture  [44],  [23]  on  the  existence  of  stable  planar  least-squares  inverses  for 
a  given  2-D  polynomial,  and  on  the  existence  of  2-D  Beurling  factorizations  in 
the  Hardy  space  H^  (cf.  [45]). 


VI .  Conclusions  and  Extensions 

In  this  paper,  we  have  presented  some  recursive  algorithms  for  the  deter¬ 
ministic  and  stochastic  modeling  of  2-D  systems.  These  results  can  be  extended 
in  several  directions.  One  such  extension  would  be  to  consider  the  matrix  form 
of  the  results  presented  here.  We  note  in  that  respect  that  a  matrix  version 
of  the  2-D  modified  approximation  problem  of  Section  III  has  recently  been 
studied  by  Bose  and  Basu  [46].  Also,  as  noted  in  Sections  III  and  V,  the  2-D 
Lanczos  and  Levinson  recursions  that  we  have  obtained  for  the  solution  of 
Toeplitz  block  Toeplitz  equations  require  0(n  m  )  or  0(n  m  )  operations  to  invert 

a  matrix  with  m+1  blocks  of  size  n+1  x  n+1.  By  comparison,  the  1-D  Levinson  al- 

2 

gorithm  requires  only  0(n  )  operations  to  invert  an  n+1  x  n+1  Toeplitz  matrix, 

so  that  the  2-D  recursions  do  not  preserve  the  efficiency  of  the  1-D  algorithm 

2  2 

in  both  dimensions.  To  obtain  algorithms  of  complexity  0(n  m  )  (or  less),  it  is 

likely  that  the  recent  doubling  algorithms  of  Morf  [33]  and  Gustavson  and  Yun 

[32]  will  play  a  significant  role.  Another  direction  of  generalization  would 

be  to  consider  2-D  extrapolation  problems  for  the  case  when  the  process  y(. ,  .) 

is  not  only  translation  invariant,  but  also  isotropic,  i.e.,  when  the  covariance 

2  2  1/2 

of  y(t,  s)  and  y(t',  s')  depends  only  on  d  =  ((t-f)  +  (s-s')  )  ,  the  distance 

between  (t,s)  and  (t',s').  For  such  processes,  some  extrapolation  problems  pre¬ 
senting  a  circular  symmetry  have  been  considered  by  Popov  [47]  and  Yadrenko 
[48].  However,  the  extrapolation  algorithms  obtained  in  this  context  are  still 
quite  inefficient,  and  the  problem  of  finding  some  efficient  algorithms  seems 
worth  considering. 
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Figure  Captions 


Figure  1. 
Figure  2. 
Figure  3. 


Figure  4. 


A  geometry  of  modified  approximation. 

Some  geometries  of  rational  approximation. 

(a)  The  domain  of  approximation  of  k1  /h1 

n>m  n,m 

(b)  The  domain  of  approximation  of  u^  /v-5 

n  j  id  n  $  in 

(a)  The  prediction  geometry  for  the  horizontal  predictor  h* 

(b)  The  prediction  geometry  for  the  vertical  predictor  m> 


